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Abstract 



We present results of a large-scale simulation for the Kaon B parameter Bk 
in quenched lattice QCD with the Kogut-Susskind quark action. Calculating 
Bk at 1 % statistical accuracy for seven values of lattice spacing in the range 
a as 0.24 — 0.04 fm on lattices up to 56 3 x 96, we verify a quadratic a depen- 
dence of Bk theoretically predicted. Strong indications are found that, with 
our level of accuracy, ajjg{l/a) 2 terms arising from our one- loop matching 
procedure have to be included in the continuum extrapolation. We present 
-Bft-(NDR, 2 GeV)=0. 628(42) as our final value, as obtained by a fit including 
the otjfg(l/a) 2 term. 
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The knowledge of the kaon B parameter Bk 
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is imperative to extract the CP violation parameter of the Cabibbo-Kobayashi-Maskawa 
matrix from experiment. Work has been continued for a decade to determine this parameter 
with lattice QCD |1[ using both Wilson and Kogut-Susskind quark actions. Calculations 
with the latter have the advantage that the correct chiral behavior of the matrix element 
is ensured by U(l) chiral symmetry. Nonetheless, previous studies with this action 
have not yielded a definitive result for the matrix element. 

A major difficulty, uncovered in Ref. @, is the presence of a large scaling violation in 
Bk, which renders a reliable extrapolation to the continuum limit nontrivial. Whereas the 
scaling violation is theoretically expected to be 0(a 2 ) || with the Kogut-Susskind action, 
simulations so far could not confirm it due to large statistical errors. 

Another problem concerns systematic uncertainties in the renormalization factors needed 
to match the lattice result to that in the continuum. While an earlier study [[| found that 
one-loop perturbation theory is reasonably accurate, the problem of the systematic error 
associated with renormalization has not been fully explored yet. 

In order to resolve these problems, we have carried out a large-scale simulation for Bk 
with the Kogut-Susskind quark action in quenched lattice QCD. In this article we report on 
the continuum limit of Bk, expounding the crucial points of our simulations and analysis. 

The parameters employed in our simulations are summarized in Table |I[ In order to study 
the continuum limit, seven values of the inverse gauge coupling constant (3 = 6/g 2 spanning 
the range (3 = 5.7 — 6.65 are chosen for the simulations, corresponding to the lattice spacing 
a ~ 0.24 — 0.04 fm. We set the physical scale of lattice spacing by p meson mass in the VT 
channel. The physical lattice size is kept approximately constant at La ~ 2.3 — 2.5 fm in 
order to distinguish scaling violation effects from those of finite lattice. Finite-size effects are 
examined separately at f3 = 6.0 and 6.4, varying the lattice size over the range La ~ 1.8 — 3.1 
fm. Numerical simulations have been carried out on the Fujitsu VPP500/80 supercomputer 
at KEK. 

We employ both gauge-invariant and noninvariant four-quark operators H, which differ 
by an insertion of gauge link factors connecting the quark fields spread over a 2 4 hypercube. 
The bare lattice operators are mean-field improved through a replacement x ~^ y/^oX f° r 
the quark field and — > u^U^ for the gluon field, where uq = P 1 ^ 4 0] P being the average 
value of the plaquette. 

The matching of Bk between lattice and continuum is made in the following way. We 
first correct lattice values of Bk by the one-loop renormalization factor ||^] evaluated with 
the MS coupling cejjg(q*) at a matching scale q* = 1/a [0)011 to obtain the continuum 
operator S^(NDR, q*) renormalized in the MS scheme with the naive dimensional regular- 
ization(NDR). The continuum value at a physical scale \l — 2 GeV is then obtained via a 
2-loop running of the continuum renormalization group starting from 5^(NDR, q*): 
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where f3 = 11, Pi = 102, 7o = 4 and ji = —7 [12| are the Nj = quenched values for the 
renormalization group coefficients. This procedure leaves an uncertainty of 0{ajfg{q*) 2 ) in 
i?(NDR, q*) arising from the use of one-loop renormalization factors fl5,13"l. 

The coupling constant otjjg{q*) needed in the matching factor is obtained once the A-^ is 
specified. To estimate this, we start from ap |L4| defined by — logP = 47r/3ap(3.40/a)(l — 
1.19ap), and calculate Ajjg = 0.625Ap from ap(3.40/a), where the 3-loop correction term 
is included. The value of estimated in this way, however, suffers from scaling violation. 
We therefore extrapolate the results at our seven values of (3 quadratically in m p a to the 
continuum limit, finding A^- = 232(4) MeV. We then take Ajj$ = 230 MeV, and calculate 
the MS running coupling to 3-loop accuracy, which is used throughout our analyses to 
minimize additional scaling violation entering into the Bk calculation. 

In our simulations gauge configurations are generated with the 5-hit heatbath algorithm, 
and B K is calculated at every 1000(0 = 5.7), 2000(/3 < 6.0) or 5000(/3 > 6.2) sweep intervals. 
Our main results are based on calculations at four values of degenerate strange and down 
quark mass m q a, equally spaced in the interval given in Table |. 

Lattice values of Bk are calculated from the 3-point Green function of the four-quark 
operator at time t with two kaons created at the temporal edges of the lattice, divided by the 
vacuum saturation of the same operator. Eight wall sources corresponding to the corners 
of a spatial cube are employed to construct a quark-antiquark propagator combination such 
that only the pseudoscalar meson in the Nambu-Goldstone channel propagates jnj. Quark 



propagators are calculated with the Dirichlet boundary condition in time and the periodic 
boundary condition in space. Gauge configurations are fixed to the Landau gauge. 

The fitting interval to extract Bk from the Green function is chosen so that the minimum 
time t m i n a is approximately constant at t m i n a ~ 1.4 — 1.5 fm for all values of f3. The resulting 
Bk changes less than ±0.3 % for all m q a and (3 under a variation of t min by ±2. 

At each value of (3 lattice results are interpolated in m q a with the formula suggested by 
chiral perturbation theory [f[6|, 



B K — B[l — 3c • m q a\og{m q a) + b ■ m q a\. (3) 

The physical value of Bk is obtained at half the strange quark mass m s a/2, estimated from 
experiment, m K /m p = 0.498/0.770. 

We present our results in Table [TI|. The errors are estimated by a single elimination 
jackknife procedure. Our statistical error is small, being 0.1% at /3 = 5.7 and gradually 
increasing to 1.2% at (3 — 6.65. At (3 — 6.0 and 6.4, three spatial sizes are examined for a 
finite-size study. Some size dependence of order 2% is seen below the spatial size La «2.0 fm 
at (3 = 6.4, but the magnitude decreases to less than 0.5% for fm at both values of 

(3. We have made our main runs with a spatial size larger than La ~ 2.3 fm, thus expecting 
finite-lattice corrections being smaller than the statistical error. 

We present _Br-(NDR, 2GeV) as a function of m p a in Fig. [I] both for gauge noninvariant 
(circles) and invariant (diamonds) operators. The five points below m p a ~ 0.6(/3 > 5.93) are 
consistent with the 0(a 2 ) scaling behavior theoretically expected 0. Toward large lattice 
spacings, however, we observe a change of curvature from a positive to a negative sign. At 
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an intermediate range m p a ~ 0.6 — 0.3 {(3 = 5.85 — 6.2) a cancellation among the a 2 and 
higher order terms conspire to yield an apparently linear dependence of Bk- This is the 
linear behavior we observed at an early stage of our work ||17|| . The later result at a smaller 
lattice spacing m p a ~ 0.22(/3 = 6.4) gave a first indication of an 0(a 2 ) behavior |18|]; this 
is now confirmed by the calculation at a yet smaller lattice spacing m p a « 0.16(/3 = 6.65) 
given in this paper. 

In our preliminary report (TTSJ we took a naive approach to estimate the continuum Bk, 
simply by applying a polynomial fit assuming 0(a 2 ) dependence. A fit of the five points 
above (3 = 5.93 with the form B K = c + Ci(m p a) 2 , shown by the dashed lines in Fig. [TJ, gives 
a value at the continuum i?x(NDR, 2 GeV)=0. 616(5) for the gauge noninvariant operator, 
and 0.580(5) for the invariant one, the average of the two being 0.598(5). 

An obvious problem with this analysis is that the two operators yield different values. 
We recall that Bk for the two operators, and hence also their difference, should receive not 
only 0(a 2 ) scaling violation but also Oi-jfg(q*) 2 errors from the matching procedure. Figure^ 
plots the difference as a function of m p a (numerical values given in Table [n]). Small errors 
of 3-4% result from a strong correlation between the matrix elements of the two operators. 
We find that the difference can be fitted by the form bi{m p a) 2 + b 2 ajjg(q*) 2 : employing five 
data points for m p a^0.5 we obtain b\ = —0.23(2) and b 2 = 1.73(5) for x 2 /d.o.f = 2.2. The 
solid line indicates the fit, and the others show the breakdown into the a 2 (dotted line) and 
a 2 (dashed line) contributions. A fit allowing a constant b yields a value of bo vanishing 
within 2cr: bo = —0.032(16), b\ = —0.44(11) and b 2 = 3.4(8). These results strongly indicate 
that the decrease of the difference toward small lattice spacings seen in Fig. is actually an 
"msO*) 2 effect. 

Encouraged by this analysis we attempt to fit the five points at (3 > 5.93 simultaneously 
for both operators including their correlations with the form B™ n ~ mv = Co+Cia 2 +c 2 ajjg(q*) 2 
and B™ = d + d^a 2 + d 2 ajj^(q*) 2 . This yields c = 0.67(6) and d = 0.71(7), and hence 
we impose the constraint Co = do in our final fit. In the continuum limit the fit (solid lines 
in Fig. D gives B K (NDR, 2 GeV) = 0.628(42) with x 2 /d.o.f=1.37. The error is roughly 
ten times the one from the naive quadratic fit. This large error reflects uncertainties of the 
coefficient of the a 2 terms: c 2 = —0.5(2.0) and d 2 = —2.2(2.0). The difference, however, is 
well constrained; c 2 — d 2 = 1.7 agrees well with b 2 = 1.73 obtained above. 

We find larger coefficients c 2 = —1.0(4.2) and d 2 = —4.3(4.2) when q* = ir/a is used, or 
c 2 = 1.6(1.5) and d 2 = —3.2(1.5) if mean-field improvement is not made for the operators. 
This supports the tadpole argument of Ref. [/J. 

The final value depends only weakly on our choice of A-^ = 230 MeV: e.g., i?^(NDR, 
2 GeV)=0.627(42) for A I7 ^=220 MeV and 0.628(41) for 240 MeV. 

As our final value of Bk in the continuum limit, we adopt the fit including the a 2 term, 

5 X (NDR, 2GeV) = 0.628 ± 0.042, (4) 

which includes a systematic error from the 2-loop uncertainty. The size of the quoted 
error is 6.6%, which roughly equals 3 x aj^(q* = 1 /a) 2 at our smallest lattice spacing 
1/a = 4.87GeV at f3 = 6.65 where aj^-(4.87GeV) = 0.147. This magnitude of error is 
unavoidable, even with 1% statistical accuracy at each f3 achieved in our simulation, unless 
a two-loop calculation is carried out for the lattice renormalization. 

Our final result is consistent with the JLQCD value obtained using the Wilson quark 



4 



action, 5^(NDR, 2GeV) = 0.562 ± 0.064, in which the operator mixing problem is solved 



non-perturbatively with the aid of chiral Ward identities [JT9 . 

One of the systematic errors not taken into account in our final result (HI) is the effect of 
non-degenerate strange and down quark masses m s ^ m^. Analyzing this problem is difficult 
within quenched QCD since the chiral limit rrid — > with m s ^ is expected to diverge 
due to a quenched chiral logarithm |HJ. Our attempt at a verification of the logarithmic 
divergence is also inconclusive: our results for non-degenerate quarks can be fitted quite 
well either with or without the singular term. At this stage we are not able to quote the 
magnitude of error due to the use of degenerate quark masses. 

Finally our quoted error does not include effects of sea quarks. Preliminary attempts 
suggest that the quenching error may not exceed 5% or so |H|. More extensive efforts, 
however, are clearly needed to estimate dynamical quark contributions to the Bk parameter. 
Full QCD simulations should also enable us to answer the issues with non-degenerate quark 
masses. Carrying out such calculations represent the final step toward the first-principle 
determination of the kaon B parameter. 

This work is supported by the Supercomputer Project (No. 97-15) of High Energy Accel- 
erator Research Organization (KEK), and also in part by the Grants- in- Aid of the Ministry 
of Education (Nos 08640349, 08640350, 08640404, 09246206, 09304029, 09740226). 
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FIGURES 
B K (NDR, 2GeV) vs. m p a 

q =1/a, 3-loop coupling, 5 points 
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FIG. 1. Gauge non-invariant (circles) and invariant (diamonds) Sk(NDR, 2GeV) as a function 
of m p a, together with a simultaneous fit for the two operators including a 2 term (solid lines) and 
separate fits quadratic in a(dashed lines) to the five pairs of data points for (3 > 5.93. 

AB K (NDR, 2GeV) vs. m p a 

q =1/a, 3-loop coupling, 5 points 
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FIG. 2. Difference of i?x(NDR, 2GeV) between gauge noninvariant and invariant operators as 
a function of m p a. The solid line represents a fit with a 2 and a 2 terms, while the dotted (dashed) 
line is the contribution from the a 2 (a 2 ) term. 
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TABLES 



TABLE I. Run parameters 



p 


m q a 


L 3 T 


# conf. 


m p 


o _1 (GeV) La(fm) 


train t-max 


m s a/2 


P 


5.7 


0.02-0.08 


12 3 24 


150 


0.9120(7) 


0.837(6) 


2.83 


6-16 


0.0519(8) 


r\ f A r\r\r\ 

0.54900 


5.85 


0.01-0.04 


16 3 32 


60 


0.567(13) 


1.36(3) 


2.32 


10-20 


0.0201(9) 


0.57506 


5.93 


0.01-0.04 


20 3 40 


50 


0.484(10) 


1.59(3) 


2.48 


12-26 


0.0160(6) 


0.58564 


6.0 


0.01-0.04 


24 3 64 


50 


0.410(10) 


1.88(4) 


2.52 


15-47 


0.0125(5) 


0.59374 






18 3 64 


50 


0.413(12) 


1.87(6) 


1.90 


15-47 


0.0127(7) 








32 3 48 


40 


0.383(3) 


2.01(2) 


3.14 


15-31 


0.0109(2) 




6.2 


0.005-0.02 


32 3 64 


40 


0.291(10) 


2.65(9) 


2.39 


20-42 


0.00884(57) 


0.61365 


6.4 


0.005-0.02 


40 3 96 


40 


0.222(5) 


3.47(7) 


2.28 


25-69 


0.00692(29) 


0.63065 






32 3 96 


40 


0.216(7) 


3.57(11) 


1.77 


25-69 


0.00659(41) 








48 3 96 


20 


0.219(4) 


3.52(7) 


2.69 


25-69 


0.00681(26) 




6.65 


0.004-0.016 


56 3 96 


40 


0.158(4) 


4.87(12) 


2.27 


36-58 


0.00511(25) 


0.64912 



TABLE II. Results for B^(NDR, 2GeV) at each (3 calculated with the matching scale q* = 1/a. 



p 


L 2 T 


^(NDR, 2GeV) 
non-invariant invariant 


AB K 


5.7 


12 3 24 


0.8464(7) 


0.8224(7) 


0.0240(3) 


5.85 


16 3 32 


0.7798(25) 


0.7562(25) 


0.0236(11) 


5.93 


20 3 40 


0.7522(23) 


0.7229(22) 


0.0292(8) 


6.0 


24 3 64 


0.7154(23) 


0.6826(24) 


0.0328(5) 




18 3 64 


0.7174(68) 


0.6787(68) 


0.0388(12) 




32 3 48 


0.7128(14) 


0.6790(16) 


0.0339(8) 


6.2 


32 3 64 


0.6619(48) 


0.6243(45) 


0.0376(14) 


6.4 


40 3 96 


0.6428(67) 


0.6069(69) 


0.0359(10) 




32 3 96 


0.6577(122) 


0.6126(112) 


0.0451(24) 




48 3 96 


0.6415(48) 


0.6072(51) 


0.0343(11) 


6.65 


56 3 96 


0.6350(70) 


0.6055(72) 


0.0295(10) 
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